### REPLICATION CODE FOR DUNHAM, LIEBERMAN, AND SNELL (2016)


setwd("~/Dropbox (MIT)/DLS_racerisk_PLOSms/replication_data") ##Set your own directory

library ("foreign")
library("stargazer")
library("lmtest")
library("rms")
library("ggplot2")
library("texreg")
library("car")
library("VGAM")
library("plyr")
library("psych")


######### Standard controls for regressions
controls <- "+ age + female + educ + incd + manip"

sims <- 1000 ## For bootstrapping
models <- 3

risk <- read.dta("DLS_riskdata_tot.dta") 

# create numeric emotion response items
risk$sympathy <- as.numeric(risk$sympathy)
risk$shame <- as.numeric(risk$shame)

# create total risk variable
risk$riskdiab_tot <- risk$riskdiab_you + risk$riskdiab_fam + risk$riskdiab_fri + risk$riskdiab_any
risk$riskaid_tot <- risk$riskaid_you + risk$riskaid_fam + risk$riskaid_fri + risk$riskaid_any

#recodes
risk$marstat = as.factor(risk$marstat)
risk$health12 = as.numeric(risk$health12)
risk$other12 = as.numeric(risk$other12)
risk$dishealth = as.numeric(risk$dishealth)
risk$educ <- as.numeric(risk$educ)
risk$incd <- as.numeric(risk$incd)
risk$conf_gendn <- as.numeric(risk$conf_gend2)
risk$conf_racen <- as.numeric(risk$conf_race2)
risk$conf_agen <- as.numeric(risk$conf_age2)
risk$conf_racech <- as.numeric(risk$conf_race_ch)
risk$conf_agech <- as.numeric(risk$conf_age_ch)
risk$diab_charn <- as.numeric(risk$diab_char)
risk$diab_prevn <- as.numeric(risk$diab_prev)
risk$diab_testn <- as.numeric(risk$diab_test)
risk$aids_charn <- as.numeric(risk$aids_char)
risk$aids_prevn <- as.numeric(risk$aids_prev)
risk$aids_testn <- as.numeric(risk$aids_test)


#### subset data for white and black respondents
risk_w <- subset (risk, black=="Not checked" & manip =="Pass")
risk_b <- subset (risk, black=="Checked" & manip == "Pass")



### ANALYSES AND OUTPUT
# Table 4
### Start by looking at pre-test data / How does "concern" compare at baseline?

dvarnames <- Cs(conc_asth, conc_canc, conc_deaf, conc_diab, conc_hiv, conc_infl, conc_obes, conc_blin)
dvlabs <- Cs(Asthma, Cancer, Deafness, Diabetes, HIV, Influenza, Obesity, Blindness)

modelfits.conc <- vector(length(dvarnames), mode = "list")
covars <- c("Black", "Female", "Age", "Education", "Income", "Know pers w AIDS", "Know pers w Diab", 
            "Suburban", "Rural", "Married", "Own health past yr", "Others' health past year", "Discuss health", 
            "Gay or bisexual")
names(modelfits.conc) <- dvarnames

for(i in dvarnames) {
  
  modelformula.conc <- paste(i, " ~ black + female + age + educ + incd + knowaids + knowdiab + restype + married + health12 + other12 + dishealth + gaybisex")
  modelfits.conc[[i]] <- lm(as.formula(modelformula.conc),  data = risk)
  
}


stargazer(modelfits.conc[[dvarnames[1]]], modelfits.conc[[dvarnames[2]]], modelfits.conc[[dvarnames[3]]],
          modelfits.conc[[dvarnames[4]]], modelfits.conc[[dvarnames[5]]], modelfits.conc[[dvarnames[6]]],
          modelfits.conc[[dvarnames[7]]], modelfits.conc[[dvarnames[8]]],
          style = "ajps",
          covariate.labels = covars,
          dep.var.labels.include = FALSE,
          column.labels = dvlabs,
          title = "OLS Estimates of Pre-Treatment Concern / Risk Perception",
          font.size = "tiny",
          omit.stat = c("f", "adj.rsq", "ser"),
          label = "corr",
          notes = "Each column lists the outcome variable (Concern about respective disease for American health care. Standard errors in parentheses."
)


#Tables 5 and 6
### REGRESSIONS OF THOSE IN AIDS ARM, analyzed by race group


dvarnames <- Cs(riskaid_tot, perc_aids, budg_aids, ins_aids, clicks, conf_racech, sympathy, shame)
dvlabs <- c("AIDS Risk", "AIDS budg shr", "AIDS budg inc", "AIDS insure", "Web clicks","Conf race data (chg)", "Sympathy", "Shame")

covars <- c("Treatment: Race only", "Treatment: Stigma only", "Treatment: Race + Stigma", "Disease concern" , "Age", "Female",
            "Education","Income",
            "Pass manip test")

modelfits.conc <- vector(length(dvarnames), mode = "list")

names(modelfits.conc) <- dvarnames

subs <- c("Checked" , "Not checked")
titles <- paste("OLS estimates of Treatment Effects, AIDS Condition", c("(Blacks Only)", "(Whites Only)" ))
labs <- c("blaids","whaids")

for(j in (1:length(subs))){
for(i in dvarnames) {
  
  modelformula.conc <- paste(i, " ~ as.factor(aids) + conc_hiv ", controls)
  modelfits.conc[[i]] <- lm(as.formula(modelformula.conc), data = subset(risk, black==subs[j] ))
  
}

stargazer(modelfits.conc[[dvarnames[1]]], modelfits.conc[[dvarnames[2]]], modelfits.conc[[dvarnames[3]]],
          modelfits.conc[[dvarnames[4]]], modelfits.conc[[dvarnames[5]]], modelfits.conc[[dvarnames[6]]],
          modelfits.conc[[dvarnames[7]]], modelfits.conc[[dvarnames[8]]],
          style = "ajps",
          covariate.labels = covars,
          dep.var.labels.include = FALSE,
          column.labels = dvlabs,
          title = titles[j],
          font.size = "tiny",
          omit.stat = c("f", "adj.rsq", "ser"),
          label = labs[j],
          notes = "Each column lists the outcome variable. Treatment arms are dummy variables and control group is omitted category. Standard errors in parentheses."
          
)

}



#### Tables 7 and 8
### REGRESSIONS OF THOSE IN DIABETES ARM, analyzed by race group

dvarnames <- Cs(riskdiab_tot, perc_diab, budg_diab, ins_diab, clicks, conf_racech, sympathy, shame)
dvlabs <- c("Diab Risk", "Diab budg shr", "Diab budg inc", "Diab insure", "Web clicks","Conf race data  (chg)", "Sympathy", "Shame")

modelfits.conc <- vector(length(dvarnames), mode = "list")

names(modelfits.conc) <- dvarnames

subs <- c("Checked" , "Not checked")
titles <- paste("OLS estimates of Treatment Effects, Diabetes Condition", c("(Blacks Only)", "(Whites Only)" ))
labs <- c("bldiab","whdiab")


for(j in (1:length(subs))){
  print (j)
  for(i in dvarnames) {
    
    modelformula.conc <- paste(i, " ~ as.factor(diabetes) + conc_diab ", controls)
    
        modelfits.conc[[i]] <- lm(as.formula(modelformula.conc), data = subset(risk, black==subs[j] ))
    
  }
  
  stargazer(modelfits.conc[[dvarnames[1]]], modelfits.conc[[dvarnames[2]]], modelfits.conc[[dvarnames[3]]],
            modelfits.conc[[dvarnames[4]]], modelfits.conc[[dvarnames[5]]], modelfits.conc[[dvarnames[6]]],
            modelfits.conc[[dvarnames[7]]], modelfits.conc[[dvarnames[8]]],
            style = "ajps",
           covariate.labels = covars,
            dep.var.labels.include = FALSE,
            column.labels = dvlabs,
            title = titles[j],
            font.size = "tiny",
            omit.stat = c("f", "adj.rsq", "ser"),
           label = labs[j],
           notes = "Each column lists the outcome variable. Treatment arms are dummy variables and control group is omitted category. Standard errors in parentheses."
                      
  )
  

}


##### ALL TOGETHER ####

############## BOOSTRAPPED ESTIMATES  ###################
models <- 3
pdf(file="plosfig.pdf", width = 10, height = 8)


par(mfrow = c(4,4), cex = .6)
dvarnames <- Cs(risktot, percent, budget, insurance, clicks, conf_racech, sympathy, shame)
dvarlabels <- c("Risk perception", "Percent expenditure", "Budgetary increase", "Insurance premium", 
                "Web clicks for more info", "Confidence govt data", "Sympathy", "Shame")

names(dvarlabels) <- dvarnames


for(dv in dvarnames) {
  
  ############################################################
  
  
  ## Create matrix to store output - a function of model #, sim #
  output.rs1 <- matrix(NA, nrow=sims , ncol = models)
  output.wrs1 <- matrix(NA, nrow=sims , ncol = models)
  output.drs1 <- matrix(NA, nrow=sims , ncol = models)
  
  ### BOOTSTRAP, #sim draws from dataset
  for(i in 1:sims){
    ## Sampling from data set
    rows <- sample(1:nrow(risk), nrow(risk), replace=T)
    cd.temp.risk <- risk[rows,]
    
    modelformula <- paste(dv, " ~ as.factor(treat) + black + as.factor(treat):black + female + aidsexp + concerned + as.factor(manip)")
    temp.m1 <- lm(as.formula(modelformula), data = cd.temp.risk)
    
    tempdv <- paste0 ("risk$",dv)
    tempdv2 <- eval(parse(text = tempdv))
    std.dev <- sd(tempdv2, na.rm = TRUE)
    
    
    
    bt1 <- predict(temp.m1, data.frame(treat = 4, 
                                       black = "Checked",
                                       female = 1, 
                                       aidsexp = .5, 
                                       concerned = 50,
                                       manip = "Pass",
                                       type="response"))
    
    bc1 <- predict(temp.m1, data.frame(treat = 3, 
                                       black = "Checked",
                                       female = 1, 
                                       aidsexp = .5, 
                                       concerned = 50,
                                       manip = "Pass",
                                       type="response"))
    
    bt2 <- predict(temp.m1, data.frame(treat = 2, 
                                       black = "Checked",
                                       female = 1, 
                                       aidsexp = .5, 
                                       concerned = 50,
                                       manip = "Pass",
                                       type="response"))
    
    bc2 <- predict(temp.m1, data.frame(treat = 1, 
                                       black = "Checked",
                                       female = 1, 
                                       aidsexp = .5, 
                                       concerned = 50,
                                       manip = "Pass",
                                       type="response"))
    
    wt1 <- predict(temp.m1, data.frame(treat = 4, 
                                       black = "Not checked",
                                       female = 1, 
                                       aidsexp = .5, 
                                       concerned = 50,
                                       manip = "Pass",
                                       type="response"))
    
    wc1 <- predict(temp.m1, data.frame(treat = 3, 
                                       black = "Not checked",
                                       female = 1, 
                                       aidsexp = .5, 
                                       concerned = 50,
                                       manip = "Pass",
                                       type="response"))
    
    wt2 <- predict(temp.m1, data.frame(treat= 2, 
                                       black = "Not checked",
                                       female = 1, 
                                       aidsexp = .5, 
                                       concerned = 50,
                                       manip = "Pass",
                                       type="response"))
    
    wc2 <- predict(temp.m1, data.frame(treat = 1, 
                                       black = "Not checked",
                                       female = 1, 
                                       aidsexp = .5, 
                                       concerned = 50,
                                       manip = "Pass",
                                       type="response"))
    
    output.rs1[i,1] <-  (bt2 - bc2) / std.dev #race only
    output.rs1[i,2] <-  (bc1 - bc2) / std.dev # stigma only
    output.rs1[i,3] <-  (bt1 - bc2) / std.dev # race and stigma

    output.wrs1[i,1] <-  (wt2 - wc2) / std.dev #race only
    output.wrs1[i,2] <-  (wc1 - wc2) / std.dev # stigma only
    output.wrs1[i,3] <-  (wt1 - wc2) / std.dev # race and stigma
    
    output.drs1[i,1] <-  ((bt2 - bc2) - (wt2 - wc2)) / std.dev #race only
    output.drs1[i,2] <-  ( (bc1 - bc2) - (wc1 - wc2)) / std.dev # stigma only
    output.drs1[i,3] <-  ((bt1 - bc2) -(wt1 - wc2)) / std.dev # race and stigma

    print(i)
  }
 
  ###
  
  cis.rs <- matrix(NA, nrow = models, ncol = 3)
  cis.wrs <- matrix(NA, nrow = models, ncol = 3)
  cis.drs <- matrix(NA, nrow = models, ncol = 3)
  
  cis.rs90 <- matrix(NA, nrow = models, ncol = 3)
  cis.wrs90 <- matrix(NA, nrow = models, ncol = 3)
  cis.drs90 <- matrix(NA, nrow = models, ncol = 3)
  
  for(i in 1:models){
    cis.rs[i, 1] <- mean(output.rs1[,i])
    cis.rs[i, 2] <- quantile(output.rs1[,i], 0.025)
    cis.rs[i, 3] <- quantile(output.rs1[,i], 0.975)
    
    cis.wrs[i, 1] <- mean(output.wrs1[,i])
    cis.wrs[i, 2] <- quantile(output.wrs1[,i], 0.025)
    cis.wrs[i, 3] <- quantile(output.wrs1[,i], 0.975)
    
    cis.drs[i, 1] <- mean(output.drs1[,i])
    cis.drs[i, 2] <- quantile(output.drs1[,i], 0.025)
    cis.drs[i, 3] <- quantile(output.drs1[,i], 0.975)
    
    cis.rs90[i, 2] <- quantile(output.rs1[,i], 0.05)
    cis.rs90[i, 3] <- quantile(output.rs1[,i], 0.95)
    
    cis.wrs90[i, 2] <- quantile(output.wrs1[,i], 0.05)
    cis.wrs90[i, 3] <- quantile(output.wrs1[,i], 0.95)
    
    cis.drs90[i, 2] <- quantile(output.drs1[,i], 0.05)
    cis.drs90[i, 3] <- quantile(output.drs1[,i], 0.95)
  }
  
  ##aa respondents plot
  plot(cis.rs[,1],
       ylim = c( ((-0.2 * abs(min(cis.rs[,2]))) +  min(cis.rs[,2])), ((0.2 * max(cis.rs[,3])) + max(cis.rs[,3]))),     
       xlim = c(0.5, models + 0.5),
       xlab = "",
       ylab = "",
       xaxt = "n",
       pch = 16,
       main =paste(dvarlabels[dv],"\n(African American Respondents)" ))
  
  for(i in 1:models){
    lines(c(i, i), c(cis.rs[i, 2], cis.rs[i,3]))
    lines(c((i-0.03), (i + 0.03)), c(cis.rs90[i, 2],cis.rs90[i, 2]) )
    lines(c((i-0.03), (i + 0.03)), c(cis.rs90[i, 3],cis.rs90[i, 3]) )
    
    }
  
  
  axis(tick = FALSE, 1, 1:models, c("Race\nframe", "Stigma\nframe", "Race +\nStigma"), cex.axis=1.0)
  
  
  lines(c(0.5, models + 0.5), c(0,0), lty = 2)
  
  
  ##wh respondents plot
  plot(cis.wrs[,1],
       ylim = c( ((-0.2 * abs(min(cis.wrs[,2]))) +  min(cis.wrs[,2])), ((0.2 * max(cis.wrs[,3])) + max(cis.wrs[,3]))),     
       xlim = c(0.5, models + 0.5),
       xlab = "",
       ylab = "",
       xaxt = "n",
       pch = 16,
       main =paste(dvarlabels[dv],"\n(White Respondents)" ))
  
  for(i in 1:models){
    lines(c(i, i), c(cis.wrs[i, 2], cis.wrs[i,3]))
    lines(c((i-0.03), (i + 0.03)), c(cis.wrs90[i, 2],cis.wrs90[i, 2]) )
    lines(c((i-0.03), (i + 0.03)), c(cis.wrs90[i, 3],cis.wrs90[i, 3]) )
  }
  
  
  axis(tick = FALSE, 1, 1:models, c("Race\nframe", "Stigma\nframe", "Race +\nStigma"), cex.axis=1.0)
  
  
  lines(c(0.5, models + 0.5), c(0,0), lty = 2)
  
 
}
dev.off()

##### END OF MAIN CODE ##################





